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Introduction 


When solutions of the steady Euler equations are 
sought as solutions of a stationary problem, rather than 
in the limit of some temporal evolution, the solutions 
found may not be unique. However, not every solution 
of the steady equations can be expected to represent 
a physically meaningful flow field. A physically mean- 
ingful flow field is one that can be observed in nature. 
A physically meaningful steady flow field, in addition 
to satisfying the steady equations, must also be stable. 
A stable solution, in the strictest sense, represents a 
system which, when subjected to small finite pertur- 
bations, damps them out and returns to its original 
state. In mathematical terms, such a solution is known 
as asymptotically stable in the sense of Liapunov. 

An example of a flow field with multiple solutions is 
shown in figure 1. The figure shows the Mach number 
distribution in a nozzle with an exit boundary condition 
such that a supersonic flow is established downstream 
of the sonic throat located at a: = 0 in the figure. There 
are two solutions consistent with the steady equations 
(quasi-one-dimensional Euler equations in this case) 
and the given boundary conditions: one with a shock 
wave in the divergent section of the nozzle (curve abode) 
and the other with a shock wave in the convergent sec- 
tion (curve abfde). It is well-known (see ref. 1) that the 
solution with the shock in the divergent section is sta- 
ble, while the solution with the shock in the convergent 
section is unstable. The arguments given in most text- 
books (see ref. 1) to explain the instability are, however, 
not rigorous. The purpose of this report is to present a 
simple mathematical procedure to analyze the stability 
of flows involving planar shock waves. A planar shock 
wave here means one whose line of interception with 
the plane on which the solution is sought is a straight 
line. The mathematical procedure is based on condi- 
tions only in the immediate neighborhood of the shock, 
so it is therefore referred to as a “local analysis.” The 
local analysis provides necessary conditions for stabil- 
ity. To establish necessary and sufficient conditions for 
stability, a global analysis, which takes far-field bound- 
ary conditions into account, is required. The disadvan- 
tage of the global analysis is that, in general, it requires 
numerical integration to solve the resulting eigenvalue 
problem. On the other hand, with the local analysis, a 
closed-form solution can be obtained. 

The stability of a shock wave in the divergent- 
convergent section of a nozzle was studied analytically 
by Kantrowitz (ref. 2) in 1947, and numerically by 
Moretti (ref. 3) in 1971. Recently, Culick and Rogers 
(ref. 4) have repeated the analysis of Kantrowitz to 
study the response of a normal shock in a ramjet engine. 
Their analysis, like that of Kantrowitz is based, on the 
assumption of quasi- steady behavior. As a justification 
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Figure 1. Nozzle flow field with two possible solutions: one 
corresponding to a shock in the divergent section (abcde) 
and the other corresponding to a shock in the convergent 
section (abfde). 

for the quasi- steady approximation, they consider the 
ratio of the time derivative to the space derivative 
occurring in the one- dimensional continuity equation, 
namely 


^ Ap Ax 

at ^ p^ 

dpu Ap I Au 

dx p ' u 


( 1 ) 


and they assume that the characteristic length Ax 
over which spacial changes are significant is of the 
order of the shock-wave thickness. Because the shock 
thickness is of the order of a few mean free paths, they 
concluded that the time derivative may be ignored. 
There are two errors in this analysis. First, the ratio in 
equation (1) obviously has a value of —1. They should 
have considered the ratio of the unsteady term to one 
of the steady terms, not to both; for example, the ratio 


It ^ 1:^ 

u At 


( 2 ) 


The second, and perhaps a worse error, is that the char- 
acteristic length is not of the order of the shock- wave 
thickness, but rather of the order of the shock-wave ex- 
cursion. Therefore, the term Ax/At in equation (2) is 



of the order of the shock speed w. The time derivatives 
are, therefore, of the same order as the shock speed and 
are retained in the present analysis. 

The analysis presented here is first developed for 
the nozzle problem just discussed. It is then used to 
investigate how the stability of the shock changes if the 
shock is assumed to be isentropic. Finally, it is used to 
study the stability of weak and strong shocks attached 
to wedges and cones — a problem previously investigated 
by the author in reference 5. 

Symbols 

a speed of sound 

A cross-sectional area of quasi- one-dimensional 

nozzle 

c damping coefficient 

F steady-state part of compatibility equation 

j equal to 0 for a two-dimensional wedge and to 
1 for an axisymmetric circular cone 

k spring coefficient 

m mass coefficient 

M Mach number 

p pressure 

P defined by equation (17), (39), or (57) 

Q defined by equation (18), (40), or (58) 

r radial coordinate (fig. 7) 

To radial coordinate at which governing equation 

is evaluated 

R defined by equation (19), (41), or (59) 

S defined by equation (20), (42), or (60) 

t time in fixed frame of reference 

T defined by equation (64) 

u longitudinal velocity component for nozzle 

problem; circumferential velocity component 
for wedge or cone problem 

V radial velocity component for wedge or cone 
problem 

Voo magnitude of free-stream velocity 

w shock-wave speed 

X longitudinal coordinate in fixed frame of 
reference 

^ shock location relative to initial position 
Oi defined by equation (8) 


P defined by equation (14) 

7 isentropic exponent 

S defined by equation (13) 

^ damf!)ing ratio 

0 circumferential coordinate (fig. 7) 

A characteristic slope 

^ space coordinate in moving frame of reference 

p density 

T time coordinate in moving frame of reference 

w undamped natural frequency 

Subscripts: 

s value at shock wave 

2 : differentiation with respect to 

1 differentiation with respect to ^ 

r differentiation with respect to r 

0 value at i = 0 

1 low-pressure side of shock wave 

2 high-pressure side of shock wave 

00 free-stream value 

Superscripts: 

+ left-running characteristic for wedge or cone 
problem 

— left-running characteristic for nozzle problem 

A dot over a symbol denotes the derivative with 
respect to time. 

Details of the Method 

Nozzle Problem With Rankine-Hugoniot Shock 

Consider the steady flow that develops in a nozzle 
as shown in figure 1. At time t = 0, assume that the 
shock is perturbed from its steady-state position by as- 
signing to the shock some “small” velocity wq. The 
ensuing evolution will be studied to determine whether 
the shock will return to its steady-state position. With 
the assumption that the flow is governed by the quasi- 
one-dimensional Euler equations, the solution is known 
exactly upstream of the shock as a function of the nozzle 
area. Downstream of the shock, the left-running char- 
acteristic brings information from the subsonic region 
to the shock, as indicated in figure 2. The shock mo- 
tion is governed by the Rankine-Hugoniot jumps and by 
the compatibility relation reaching the shock along the 
left-running characteristic. In terms of space and time 
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Figure 2. Characteristic diagram in the neighborhood of a 
shock wave. 

coordinates (^,t) in a frame of reference that moves 
with the shock, 

r = i 1 

^ = X - Xg{t) j 

the compatibility relation valid on the left- running char- 
acteristic is 



P 2 

Pi 

( 10 ) 

U2-W l + ^(mi-^)2 

Ui-w 

( 11 ) 

where 


1 ( 
II 

( 12 ) 

2 

(13) 


(14) 

Equations (10) and (11) are valid along the shock path, 
^ = 0. Therefore, by taking the derivative of 

these two equations, the rate of change in time of p 2 
and U 2 can be related to the changes taking place on the 
supersonic side of the shock due to the shock motion. 
After performing this operation, the resulting equations 
are linearized by neglecting terms of order {w/ax)^ 
and higher. Finally, the derivatives dMx/dr^dpx/dr, 
and daifdr are replaced by the quasi-one-dimensional 
relations given in appendix A. The final expressions are 

— = Paw -b Qwt 

P 2 

(15) 

— = Raw + bwr 

(16) 


where 


= u — a 

(6) 

M=- 

(7) 

a 

1 dA{x) 
A{x) dx 

( 8 ) 


(9) 


In the above equations, (x^t) are the longitudinal and 
time coordinates in the fixed frame of reference; p, p, 
and a are the density, pressure, and speed of sound; 
u is the velocity in the longitudinal direction; 7 is the 
isentropic exponent; Xs is the shock position; and A(x) 
is the cross-sectional area of the nozzle. 

The jump conditions connecting the supersonic side, 
denoted by subscript 1 , to the subsonic side, denoted by 
subscript 2 , are 



With the aid of equations (15) and (16), equation (4) is 
evaluated at the shock to get 


(7M25 - Q)wr 

+ qM2R-P-\-^ ^ aw-F^O ( 21 ) 

M.2 — 1 
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where 


F = A H- ^ua (22) 

Equation (21) describes the shock motion. In the 
analysis of references 2 and 4, the equation describing 
the shock motion is obtained, somewhat arbitrarily, 
from one of the shock jump conditions. Equation (21), 
on the other hand, simultaneously accounts for the 
interaction of the pressure and velocity fields on the 
shock motion. In the steady state, F = 0. Considering 
departures from equilibrium to be small allows the term 
F appearing in the coefficient of w in equation (21) to 
be neglected. With the notation, 

m = 7M25 - Q (23) 


c=(pM,R-P+^y (24) 

equation (21) can be written as 


mwr -\-cw = F (25) 

The solution to this equation is given by 


w(t) — exp ( — / ~dt] 

\ Jo rn J 

[/ “‘’(i ^ 



Figure 3. Behavior of m,c/a, and k/a^ as functions of Mi. 
7 = 7 / 5 . 

where the derivatives in equations (29) and (30) are, of 
course, evaluated at the subsonic side of the shock and 
r = 0 or = 0. The first order, or linear, approximation 
to equation (25) is, therefore, 

mz cz kz = 0 (31) 

where 


If F is neglected entirely, the solution simplifies to 

w{t) = ^^;o exp J ~ (27) 

As shown subsequently, the coefficients m and c/a are 
positive for all values of Mi. Therefore, the shock 
velocity given by equation (27) grows exponentially if 
O' < 0 or decays exponentially if a > 0. This proves 
stability in the sense of Liapunov (ref. 6). 

Now, the effect of F on the motion of the shock is 
taken into consideration. Let 2? be the displacement of 
the shock wave from its equilibrium position, 

z{t) = Xs{r) - Xs^Q ~ f wdt (28) 

Since only small deviations from equilibrium are con- 
sidered, F may be expanded in a Taylor series, 

F{t) = FrT + Frr^ + ■ ■ . (29) 

or alternatively, 

F{z) = F^z + F^^^ + ■ ■ ■ (30) 


F 

k = -F, = (32) 

Wo 

and where, to simplify the notation, a dot denotes dif- 
ferentiation with respect to r. According to Liapunov’s 
theorem (ref. 6), the information obtained from the lin- 
ear approximation (i.e., eq. (31)) is sufficient to give a 
correct answer to the question of stability of the non- 
linear equation (i.e., eq. (25)) as long as F^ is nonzero. 
If the coefficients of equation (31) are assumed to be 
constant, then the equation describes the well-known 
damped harmonic oscillator; see, for example, refer- 
ence 7 or 8. The coefficients m, c, and k are the re- 
spective mass, damping, and spring coefficients of an 
equivalent oscillator. After several simplifying assump- 
tions, discussed in appendix B, the spring constant may 
be written as (eq. (BIO)) 

^^2 / c \ 

^ "" I-M2 (m) " 

The behavior of m,c/Q;, and k/a"^ as functions of the 
upstream Mach number M\ is shown in figure 3 for 
7 = 7/5. The characteristic behavior of equation (31) 
is governed by the damping ratio 
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and the undamped natural frequency, 



Since c is proportional to a, k is proportional to The 
spring coefficient k is, therefore, independent of the sign 
of a. The damping ratio on the other hand, has the 
same sign as a, but is independent of the magnitude of 
a. 

In general, the following five cases are recognized: 

(1) ^ < 0 Unstable 

(2) ^ — 0 Undamped 

(3) 0 < ^ < 1 Underdamped 

(4) ^ = 1 Critically damped 

(5) ^ > 0 Overdamped 

The behavior of ^ for a > 0 as a function of 
Ml is plotted in figure 4. The local analysis predicts 
an underdamped behavior for the shock motion in 
the divergent section. As Mi approaches infinity, f 
approaches a value of 0.74. In the convergent section 
{a < 0), the magnitude of f is the same as that 
shown in figure 4; however, its sign is negative and, 
therefore, the shock is unstable. The behavior of the 
undamped natural frequency is shown in figure 5. The 
case a = 0 corresponds to case (2), and the system is 
said to be undamped. However, for this case, the linear 
approximation does not meet the Liapunov criterion 
that Fz be nonzero. Thus, further analysis is required 
to make a conclusive statement about this case. 



Figure 4. Behavior of damping ratio ^ as a function of Mi. 



Figure 5. Behavior of undamped natural frequency as a 

function of Mi . 

Nozzle Problem With Isentropic Shock 

For weak shocks, it is generally accepted that the en- 
tropy generated at the shock can be neglected without 
seriously compromising the quality of the results. Un- 
der this approximation, the velocity field is described 
by the gradient of a potential function that satisfies 
the conservation of mass equation. The pressure and 
density are obtained by satisfying the conservation of 
energy equation and the isentropic pressure-density re- 
lation. For the nozzle problem illustrated in figure 1, 
the potential approximation to a shock connects the su- 
personic branch (curve abf) to the isentropic subsonic 
branch (curve agh). Therefore, there is an indetermi- 
nacy in the shock position, since a shock anywhere in 
the channel satisfies the imposed back pressure corre- 
sponding to point h in figure 1. The problem is dis- 
cussed in detail in reference 9. The interest here is 
to examine the stability of such an isentropic shock 
which, in view of the above remarks, is expected to be- 
have differently from the Rankine-Hugoniot shock just 
investigated. 

For this problem, the governing equations are the 
time- dependent conservation of mass and energy equa- 
tions and the isentropic pressure-density relation. With 
these equations, it is easy to show that once again equa- 
tion (4) is valid along the left- running characteristic 
reaching the subsonic side of the shock. The jump re- 
lations across the shock are now given by 
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Px{ui -w) = P2{U2 - w) 


(36) 


+ 6(ui - + 6{u 2 - w)^ (37) 

Pi P2 

- = (38) 

Pi \PiJ 

from which, proceeding as before, one gets equa- 
tions (15) and (16) but with P, Q, P, and S now defined 
as 



Now evaluating equation (4) at the shock results in an 
equation equivalent to equation (31), namely, 


where 


mz — 0 

(43) 

m = 7 M 25 - Q 

(44) 

c = 0 

(45) 

k = 0 

(46) 


For the isentropic shock, the damping coefficient is 
identically zero. The shock is, therefore, undamped. 
However, like the Rankine-Hugoniot shock with a = 0, 
the linear approximation fails the Liapunov criterion, 
and further analysis is required to study the effect of F 
on the solution. The isentropic and Rankine-Hugoniot 
mass coefficients are compared in figure 6. 


Planar Shock Attached to a Wedge or Cone 

For flow over a wedge or cone whose deflection angle 
is less than the angle associated with shock detachment 
for an incoming supersonic flow, the steady-state Euler 
equations admit two different solutions. The solutions 
are labeled according to the strength of the attached 
shock, which can be either strong or weak, depending 
upon whether its inclination is greater or less than the 
shock- wave inclination at detachment. The stability of 



Figure 6. Comparison between isentropic and Rankine- 
Hugoniot mass coefficients as a function of M \ . 



Figure 7 . Diagram of coordinate system and velocity com- 
ponents used for shock attached to a wedge or cone. 

these two solutions has been studied numerically in ref- 
erence 5, to which the reader is referred for a derivation 
of the governing equations and other pertinent infor- 
mation. For the present purpose, one need consider 
only the compatibility relation reaching the shock on 
the high-pressure side, namely, 

Pr \'^—W f p£ 

— h 7M 1 ( -h 7M— 

P U \p 

+ “ -h {v -f ucot^)yj = 0 (47) 

where v and u are the radial and circumferential velocity 
components as indicated by figure 7; also, j = 0 for the 
two-dimensional flow over a wedge, and i = 1 for the 
axisymmetric flow over a circular cone. Equation (47) 
is valid on the left-running characteristic where 

=u + a (48) 
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and where, as before, a frame of reference moving with 
the shock is used: 


( 60 ) 


T = t 

C = [^ - Os{t)]ro 

Here Os{t) defines the shock-wave inclination, and w = 
VoOs is the shock-wave speed. 

If one denotes by a subscript 1 conditions on the 
low-pressure side of the shock that depend on the shock 
location and by a subscript oo conditions on the low- 
pressure side that are constant, the free-stream Mach 
number is given by 




m = Voo sin^ 


(51) 


1 

1-^-2 

fs 

Poo Y 

U 2 

L P2 


P2 /. 


With z = — Osft, the equation governing the motion 

of the shock is equation (31), where now 


m = Q + 7 M 25 (61) 


and 


T 

c = P + '^M^R + -rx 




To \a2 ^ M 2 + 1 / 


(62) 

(63) 

(64) 


^1 = Voo cos 0 (52) 

The component of the free-stream Mach number normal 
and relative to the moving shock is given by 


Ml 


(ui — w) 

^00 


(53) 


The jump conditions relating the low-pressure side to 
the high-pressure side, denoted by subscript 2, are 


— ly = l?2 — ty 


^ (55, 

Poo ^ 

^2-w _ 6Ml + 1 

ui-w ml ^ ^ 

Proceeding as before, one gets equations (15) and (16) 
where and S are now defined as 

P — ^ M^ sin^cos^ (57) 


The behavior of f as a function of free-stream Mach 
number and shock-wave inclination is shown in figure 8 
for 7 = 7/5 and = 1. As expected, the local 
analysis predicts that both weak and strong shocks are 
stable. This analysis, as already mentioned, does not 
include the effect of far-field boundary conditions which, 
according to reference 5, must be included to determine 
the overall stability of the problem. The behavior of the 
undamped natural frequency is shown in figure 9. 

Conclusions 

A procedure has been developed to study the lo- 
cal stability of planar shock waves. In general, the 
equation governing the shock motion is equivalent to 
the equation of a damped harmonic oscillator. For a 
Rankine-Hugoniot shock in a divergent- convergent noz- 
zle, the analysis predicts a stable solution in the diver- 
gent section and an unstable solution in the convergent 
section, in agreement with well-known experimental ob- 
servations. For an isentropic shock in a nozzle, the anal- 
ysis predicts an undamped behavior. For shock waves 
attached to wedges and cones, it predicts a stable solu- 
tion for both strong and weak shocks. 


f Poo\ f Moosine\ 

^ /? V P2 j I aoc ] 



_ Poo Y 
V/^ P2 /. 


COS 6 (59) 
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Figure 8. Behavior of damping ratio ^ for cones and wedges Figure 9^ Behavior of undamped natural frequency u) for 

at Afoo — 2,3, and 4 as a function of the shock- wave cones and wedges at A/qq == 2,3, and 4 as a function of 

inclination. 7 = 7/5; ro = 1. the shock-wave inclination. 7 = 7/5; = 1. 


Appendix A 

Steady, Quasi-One-Dimensional Flow-Field 
Dependence On Area Variation 

The following expressions are derived in reference 1: 


dM _ 1 + 6M^ dA 

W ~ 1 - Af2 A 

du 1 dA 

T ^ “ 1 - M2 T 


(Al) 

(A2) 


da _ 6M^ dA 

T ~ 1 - M2 T 

dp _ M2 dA 
p ~ I- M'^~A 


(A3) 

(A4) 


dp _ 7M2 dA 
p “ 1 - M2 T 


(AS) 


At the shock, the following relation is also valid: 


A dr 


= aw 


(A6) 


9 



Appendix B 

Evaluation of Spring Constant k 


The above expression is to be evaluated at r = 0 on 
the subsonic side of the shock, denoted by subscript 2. 
From equation (Bl), one has 


Equation (22) may be written as 


F = — {p^ — pau^) + 7ua (Bl) 


Therefore, 
Fr^ 


^ ) (P€ - pan) — + (7wa)r 

r / r P 

(B2) 

From equation (4), it follows that 


{p^-apu^)^ 

_ (Pt - paur)^ + (A“ - w)^ (p^ - pau() + ('ypua)^ 
X~ ~ w 

(B3) 

The term (p^ - pau^)^ cannot be evaluated. However, 
in reference 10, it was shown for a number of prob- 
lems similar to the one considered here that along the 
shock path, the term (pr — paUr) is approximately zero. 
Therefore, the first term in the right-hand side of equa- 
tion (B3) is neglected. Substituting this equation in 
equation (B2) yields 


-{[(t).-(t) 

-(y) 


(A -w)^ 


— w 


X~ — w 


+ ( 7 UQ; 




(p^ - pau^) 

(B4) 


{p^ - pau^) |2,o = 


—^upa 


X- 


2,0 


and, therefore, 


■FV|2,0 “ 


~^uawj 


X- 


(B5) 


(B6) 


2,0 


But, from equation (25) evaluated at r — 0, one gets 


CWq 

'W^r,0 = 

m 

(B7) 

and, therefore, 


1 1 M 2 C 

iv 2,0 = , , , OiWQ 

M 2 — 1 m 

(B8) 

Since, 


p 1 _ 

2,0 

Wq 

(B9) 

the spring constant is given by 


, 1 M 2 c 

K = a 

l-Mam 

(BIO) 

The spring constant for the wedge or cone problem, 
given by equation (63), can be evaluated in the same 
manner by replacing ^ua in equation (Bl) by T and 


repeating the analysis. 
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